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PACS 87 . 23 . Cc - Population dynamics and ecological pattern formation 
PACS 02 . 50 . Ey - Stochastic processes 

PACS 05. 40. -a - Fluctuation phenomena, random processes, noise, and Brownian motion 
PACS 87 . 10 . Mn - Stochastic modeling 

Abstract. - Generalizing the cyclically competing three-species model (often referred to as the 
rock-paper-scissors game), we consider a simple system of population dynamics without spatial 
structures that involves four species. Unlike the previous model, the four form alliance pairs 
which resemble partnership in the game of Bridge. In a finite system with discrete stochastic 
dynamics, all but 4 of the absorbing states consist of coexistence of a partner-pair. From a master 
equation, we derive a set of mean field equations of evolution. This approach predicts complex time 
dependence of the system and that the surviving partner-pair is the one with the larger product 
of their strengths (rates of consumption). Simulations typically confirm these scenarios. Beyond 
that, much richer behavior is revealed, including complicated extinction probabilities and non- 
trivial distributions of the population ratio in the surviving pair. These discoveries naturally raise 
a number of intriguing questions, which in turn suggests a variety of future avenues of research, 
especially for more realistic models of multispecies competition in nature. 



Introduction. — Over the years evolutionary game 
theory and population dynamics have yielded important 
insights into biodiversity and the behavior of multispecies 
ecological systems [IH3] . The complexity of real- world sys- 
tems makes a full understanding of their properties very 
difficult. For that reason, the study of simple model sys- 
tems is extremely valuable, as the complete knowledge of 
these systems allows to identify generic features valid in 
the more realistic but also more complex situations. 

In this context multispecies models with cyclic compe- 
tition constitute some of the simplest cases where coex- 
istence and species extinction can be studied using tech- 
niques from statistical mechanics and from non-linear dy- 
namics [3l|4]. Many recent investigations revealed a rich 
and complex behavior. In particular, for systems with 
three species (a.k.a. rock-paper-scissors game) [5 -25 , the 
results range from surprising survival/extinction probabil- 
ities in models with no spatial structure to pattern forma- 
tion and mobility effects in one- and two-dimensional lat- 
tices. By contrast, far less is known for systems with more 
than three species 0. Frachebourg et al. [5j[6l[26] consid- 



ered M species in one-dimension, X m + X m+ i 2X m 
(m — 1, • • • , M; Im+i = Xl), competing with equal 
rates, k m . The steady states consist of single-species do- 
mains for M = 3,4, but are qualitatively different for 
M > 5. For systems on two-dimensions with slightly more 
complex rates, the segregation process and the forma- 
tion of defensive alliances have also been studied [2TH32] . 
In a recent paper [33] the Fokker-Planck equations for 
conserved quantities where derived for the M = 3 and 
M = 4 cases with equal competition rates. Motivated by 
real- world systems, some studies also focused on a large 
number of competing species with complicated interaction 
schemes [341135]. 

In this paper we investigate systematically the proper- 
ties of a non-spatial game involving four cyclically compet- 
ing species with arbitrary rates. Not surprisingly, the be- 
havior of the system is much richer than the three species 
case [TTl fTT]. For example, the number of absorbing states 
is not fixed at three (or four) but is 2 (TV + 1), where TV is 
the number of individuals in the system. Such a result can 
be intuitively understood, much like in the game of Bridge, 



1 In lattice models one sometimes speaks of a four-state rock- 
paper-scissors games when empty sites are allowed |15II24| . In the 



following we do not consider empty sites to form an independent 
species. 
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where the four players form partnerships. As a result, the 
composition of the end state tends to be coexistence of 
partner-species, e.g., L of X\ and N — L of X3. A mean 
field approximation is formulated and studied analytically. 
Capturing most of the complicated time evolution of the 
full, stochastic model, it describes much of the rich be- 
havior. Of course, it cannot predict extinction events. To 
explore those processes, we rely on computer simulations 
and discover very complex extinction scenarios that de- 
pend on both the rates and the initial conditions. Unlike 
the three species case, our system does not support 'the 
survival of weakest' or the 'law of stay-out' [17,36,37 . In- 
stead, the best maxim seems to be: "The prey of the prey 
of the weakest is the least likely to survive." This result is 
intuitively reasonable, since the prey of the weakest sur- 
vives easily and, in turn, causes its prey to die quickly. 
Notably, this maxim also applies to the three-species case, 
as illustrated by say, X\ being the weakest. Then, 'the 
prey of the prey of the weakest' is A3, the demise of which 
is excellent news for Xi! 

In the next section, we specify our model, discuss its ab- 
sorbing states, and provide analytic results in a mean field 
approach. Much of the system's qualitative behavior can 
be understood. However, predicting extinction probabili- 
ties is much more challenging and numerical simulations 
for exploring them are discussed in the following section. 
We end with a summary and some outlook for future re- 
search. 

Model specifications and mean field theory. — 

Our system consists of N individuals each of which is 
identified as one of four interacting species: A, B, C, and 
D. Endowing the species with cyclic competition, our dy- 
namics consists of picking a random pair and letting the 
interactions 

A + B V A A + A; B + C V -\B + B 
C + D^C + C; D + A 7 ^ D + D 

occur with probabilities p m , m = a, 6, c, d. Note that AC 
and BD pairs are non-interacting. Denoting the numbers 
of each species in our system by A" m , a configuration of 
the system, which has no spatial structure, is completely 
specified by these integers. With N = Y^m being 
a constant, our configuration space is actually a set of 
points within a regular tetrahedron [38 . Unlike the cyclic 
competition of three species, we have 2 (N + 1) absorbing 
states here. They form two fixed lines, N a + N c = N 
and N b + N d = N, and describe coexistence of the non- 
interacting pairs, A-C and B-D, respectively. Moreover, 
note that each face of the tetrahedron is also 'absorbing,' 
in the sense that transitions into the face are irreversible. 
Within each face, the problem is a special limit of the three 
species model, namely, one of the three rates being zero. 

From these dynamic rules, it is simple to write a mas- 
ter equation for P ({N m } ; £), the probability for find- 
ing the system t steps after an initital configuration 



{A^mo}. To find its solution is far less simple, how- 
ever. Instead, we will exploit a mean field approxima- 
tion for the evolution of the averages of the fractions, 
A (t) = E{ Nm } (MJN) P ({N m } ; t) , etc. Following stan- 
dard routes, we start from the master equation for P and 
consider the large N behavior to arrive at [39] 

d t A = [k a B - k d D] A- d t B = [k b C - k a A] B (1) 
d t C = [k c D - k b B] C; d t D = [k d A - k c C] D (2) 

where A (t) is simplified to A, etc. Here, t is regarded as 
a continuous variable and the "rates" k m can be related 
to the discrete time step, the p's above, and N [39]. Of 
course, the conservation law now reads A + B + C + D = 1. 
Since an overall scale can be absorbed into t, we will follow 
the normalization in the literature: k a + k b + k c + kd = 1. 
The remainder of this section will be devoted to a study 
of the evolution of A(t), B (£), etc. starting with A (0) = 
4) = N a0 /N, etc. 

Exploiting the exponential nature of typical 
growth/decay, we write the above equations as 

d t lnA = k a B-k d D; d t \nC = k c D-k b B (3) 
d t \nB = k b C -k a A; d t \nD = k d A-k c C (4) 

These clearly expose the alliance into opposing pairs AC 
and BD, an essential feature absent in the three species 
model. Borrowing the language of Bridge, we will refer to 
AC and BD as partner-pairs, as each player works in favor 
of its partner and against the opposing pair. To be quanti- 
tative, we construct appropriate linear combinations such 
as dt [k b In A + k a In C] = \D, with 

A k a k c - k b k d (5) 

being a crucial parameter. In addition to controlling how 
each species affects the growth/decay of the opposing pair, 
A generates a simple evolution 

Q(t) = Q{0)e xt (6) 

for the quantity 

j{k b +k c (jk d +k a 

® = B k c+ k d D^+kb ' ( 7 ) 

Similar to R = A kb B kc C ka in the three species system 
[17], Q is t-dependent as opposed to R being invariant. 
Furthermore, since A, B, C, D are bounded by unity, the 
indefinite decay/growth in Q can only occur when A, C or 
B, D vanish. As a result, the sign of A controls which pair 
survives. Intuitively, this prediction seems understand- 
able: The pair with the larger rate-product (k a k c or k b k d ) 
wins. 

Obviously, systems with A = are special, as Q is a 
constant of the motion. Indeed, there are two invariants, 
which can be simply A kb C ka and B kd D ka (as generaliza- 
tions of AC and BD in [3)133] , where k m = 1, Vra). Fixing 
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Fig. 1: An example of mean field evolution for A = 0: {k m } = 
(0.4,0.4,0.1,0.1). (a) Closed loop in the tetrahedron forming 
the configuration space and (b) the fractions of the different 
species as a function of time. These data were generated by 
applying a fourth order Runge-Kutta scheme with At = 10 -5 
on Eqns. (ITI21) . 

these by the initial conditions (Aq, Bo, Co, Do), we define 



natural variables: pa = (A/Aq) b , etc., which obey 



Pa Pc 



PbPd- 



(8) 



These equations define hyperbolic sheets through the 
tetrahedron and their intersection is a closed loop that 
resembles (the rim of) a saddle. Fig. [TJi shows an ex- 
ample of such an orbit, for the case k a = k^ = 0.4 and 
k c = kd = OA. Meanwhile, Fig. [TJd shows the associated 
ever-lasting oscillations in A, B,C, D. 
One important characteristic of a closed orbit is its 
extremal points. For example, let A± denote the 
largest /smallest values A assumes, given an initial point 
(Aq, Bo, Co, Dq) and a set of fc's. Then, A± are solutions 



to a generically transcendental equation: 



A~ kb/ka A± + 

CoA ± a = constant, which depends on Bo, Do, and the 
fc's. Typically, two distinct solutions exist, corresponding 
to the two extremes. At these points, the values assumed 
by B, C, D are, in general, not extremal themselves. While 



C = Co (Ao/A±) kb ^ ka at these points, B takes on the same 

/ 7 77 \l/(fca+fcd) 

value at both turning points: f k d a k~ ka BqDq 
Similarly, D is also unique. When A+ = A-, we are at a 
fixed line - formed by the intersection of the two planes: 
k a A — kbC and k a B — k^D. Unlike the lines of absorb- 
ing states (A-C and B-D), points on this line are neither 
stable nor stationary under the stochastic dynamics. Be- 
ing straight and bridging the A-C, B-D lines, this fixed 
line is enclosed by every closed orbit. In its neighborhood, 
these orbits approach circles, on which the system 'moves' 
with uj oc \/k a k c . Details supporting these remarks will be 
provided in a future publication [39] . 

For systems with A ^ 0, non-trivial fixed points cannot 
exist (as InQ — >• Xb). Thus, in a finite system, extinction 
of one of the species must occur quite rapidly. Fig. [2] 
shows two typical cases, one for each sign of A. Although 
the mean field provides good fits for short times, it pre- 
dicts neither the (average) time for the first species to die 
out nor the composition of the other three at this extinc- 
tion event. Nevertheless, we can again rely on mean field 
theory after the system 'lands' on an absorbing face (of 
the tetrahedron). For example, if D vanishes first and 
the system consists of (Ai, Bi,Ci,0) at that time, then 
mean field theory predicts the system to end at a point 
on the A-C line: (Af, 0, 1 — Af, 0). Here Af is the larger 
of the two solutions to another transcendental equation: 
A kb (1 - A f ) ka = A kb C ka . As will be shown below, these 
predictions are born out quite well in finite, stochastic sys- 
tems. 

Stochastic evolution: exact and simulation re- 
sults. — For a full stochastic process, there are limita- 
tions to a mean field approach. In particular, it fails when 
any N m is not 'large,' e.g., in small systems, or near ex- 
tinction events. If a system is 'very small,' numerically 
exact methods can be exploited to find exact solutions to 
the master equation. Indeed, for the smallest, non-trivial 
system (N = 4), simple algebra is sufficient for finding an- 
alytic expressions for all transition probabilities (from any 
state to any other) and for arbitrary rates. Unlike N = 3 
in three-species [17] however, the results are not trivially 
linear (in fc m ). Deferring details to elsewhere [39], we only 
present some general observations here. Though there is 
a finite probability that the 'weakest' is the lone survivor, 
there is also a good chance for (one or both of) its oppo- 
nent pair to survive. The clearest conclusion is: When 
the consumption rate of the weakest approaches zero, the 
survival probability of its partner vanishes. These consid- 
erations led us to a general maxim: "The prey of the prey 
of the weakest is the least likely to survive." As pointed 
out above, this maxim is consistent with "survival of the 
weakest" in 3-species models. 

For systems with larger TV's (say, > 100), Monte Carlo 
techniques are necessary to uncover interesting behaviors 
in our system. To speed up the runs, we exploit the 
Gillespie updating scheme, in which an interaction always 
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(b) 

Fig. 2: Examples of mean field evolution for A / 0: (a) 
{k m } = (0.45, 0.33, 0.14, 0.08) with A = .0366, and (b) {k m } = 
(0.35, 0.42, 0.09, 0.14) with A = -.0273. The initial fractions in 
both are {N m o/N} = (0.02,0.10,0.48,0.40) Note the two tra- 
jectories end on the A-C and B-D lines, respectively. These 
data were generated by the same scheme as in Fig. [T] 

occur at each 'step' (with appropriate relative probabili- 
ties) [101HT]. By contrast, the standard scheme is much 
slower, as many randomly chosen pairs do not interact. Of 
course, the detailed t dependence will be quite different, 
so that direct comparisons with mean field predictions are 
not possible. Nevertheless, we can rely on this scheme 
for drawing conclusions on survival and coexistence (as 
we can show that the extinction probabilities are scheme- 
independent). In particular, starting with many random 
initial conditions (typically 20000) and a variety of rates, 
simulations with N = 100 K confirm the mean field pre- 
dictions, namely, runs for A ^ systems ending with the 
correct partner pairs and runs for A = cases failing to 
end [42]. In the remainder of this letter, we will focus on 
a few particular systems with intermediate TV's, thereby 
emphasizing extinction processes and highlighting the dif- 
ferences between stochastic and mean field trajectories. 



As expected, this difference is most pronounced for 
A = systems. Instead of closed orbits (blue online in 
figures), the trajectories in a finite stochastic system (red 
online) end in an absorbing state. In the 3-species case, 
they mostly end with the 'weakest' species as sole sur- 
vivors. In our model, no simple conclusions can be drawn. 
Figs. [3^i,b illustrate a good example, where runs with iden- 
tical initial conditions end very differently. Specifically, we 
have N = IK, initial fractions (0.02, 0.10, 0.48, 0.40) and 
rates (0.4, 0.4, 0.1, 0.1). The stochastic trajectories follow 
these closed loops closely at early times. But, the noise 
drives Q away from Q (0), so that they later diverge sig- 
nificantly and, after one of the four species dies out, end 
rapidly on an absorbing state. Though both runs have the 
same {A/" m o,/c m }, the final states in Fig. [3^i and b con- 
sist of opposite partner-pairs: AC and BD, respectively. 
In addition, there are non-trivial distributions of survival 
fractions within each pair (e.g., Af, when AC survives). 
Further details will be published elsewhere [39] . 

Systems with A ^ may appear uninteresting, as they 
evolve quickly towards absorbing states. However, we 
discover rather complex behavior, especially for systems 
with extreme rates. Let us provide one illustration, with 
N = IK, initial fractions (0.1, 0.7, 0.1, 0.1), rates (0.1, 
0.0001, 0.1, 0.7999) and 10 K independent trails. Since 
A > 0, the survival rate of the weakest (B) is low (~ 10%), 
while 90% of the runs end on the AC line. Fig. [4^ 
shows one particular stochastic trajectory (red online) in 
the latter class, as well as the mean field orbit (blue on- 
line). Note that they come very close to the ABC face, 
(i.e., D <C 1) in two earlier occasions. Not surprisingly, 
in these close encounters, many runs actually 'land' on 
this face. Fig. 0}d displays the 'landing sites' (Ai,C{ with 
Bi = 1 — Ai — d) from these runs. Notably, they fall 
into three clusters (red, black, green online). Predicting 
the remarkable shape of the black cluster will undoubt- 
edly be a serious challenge! From here, the system quickly 
evolves to the A-C line, into similarly colored clusters. As- 
sociated with the unusual black cluster shape, we find the 
distribution of the final Af to be highly skewed and non- 
Gaussian. More details, as well as possible explanations, 
will be provided in a later publication. Here, let us focus 
on the evolution from a 'landing site' (Ai, Cj) to the final 
point (Af,Cf). For the ~ 9K runs in this series, we com- 
pute k = In (Af /Ai) + 1000 In (Cf/Ci) as a sensitive test of 
the mean field prediction that (A f /Ai) kb (C f /Ci) ka = 1. 
Though not identically zero, less than 4% of the values of 
k (out of the ~ 9K values) are outside the range [—1,1]! 
It is fair to conclude that, in this respect, the mean field 
approach is extremely successful. 

Summary and outlook. — In this article, we inves- 
tigate the time dependent behavior and extinction proba- 
bilities of a simple model of population dynamics: TV in- 
dividuals of four different 'species,' competing cyclically. 
Though seemingly a trivial extension from a similar three- 
species game (rock-paper-scissors), this system displays 
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Fig. 3: Two examples of stochastic evolution (thin dashed 
lines, red online) of the system in Fig. [T] both with {Nmo} = 
(0.02, 0.10, 0.48, 0.40) x 1000. They follow the mean field loop 
(thick line, blue online) initially, but diverge eventually, ending 
with opposite partner-pairs: (a) on the A-C line and (b) on the 
B-D line. 



much richer phenomena. Since the configuration space 
here is (the interior of) a tetrahedron rather than a tri- 
angle, trajectories of the system may twist and turn in 
3-d. Since the four form 'partner pairs', much like in the 
game of Bridge, the end states typically consists of one 
of the pairs, with TV — 1 possible compositions in each 
case. As a result, there are 2 (TV + 1) absorbing states 
(instead of just 3 or 4), with generally non-trivial distri- 
butions among them. The faces of the tetrahedron are 
also 'absorbing' (in that they correspond to the extinction 
of at least one species); yet the trajectories on them are 
not just trivial straight lines. Concerning extinction sce- 
narios, a law gleaned from previous studies - 'survival of 
the weakest' - seems to be violated here. Instead, our ob- 
servations support, most consistently, a different maxim: 
"The prey of the prey of the weakest is the least likely to 
survive." Much easier to understand at the intuitive level, 



B 




(a) 




(b) 



Fig. 4: A system with 'extreme' rates, {k m } — 
(0.1,0.0001,0.1,0.7999), starting with {N m0 } = 
(0.1,0.7,0.1,0.1). (a) One stochastic trajectory (thin 
dashed line, red online) and the mean field evolution (thick 
line, blue online), (b) Scatter plot of the composition in 8870 
runs of the system at the moment D becomes extinct (i.e., 
'landing sites' on the ABC face), shown as three separate 
clusters (red, black, green online) with A + C < 1. Each run 
ends on an absorbing state (a point on the A + C — 1 line). 
Also shown is the scatter plot of these associated states, with 
the three clusters labeled, respectively, by 1,2, and 3. 



this maxim also applies to the three-species game, where 
the demise of the prey of one's prey also enhances one's 
survival! 

Using a mean field approach and computer simulations, 
we report a number of other notable findings. Similar 
to, but more interesting than, R in [17], our quantity Q 
(Eqn. [7j) grows/decays exponentially in mean field theory 
and serves as an excellent indicator for which partner pair 
will survive. This theory is also quite successful in pre- 
dicting the evolution of the stochastic system, as long as 
(a) no species is close to extinction, or (b) one species is 
extinct. These two seemingly contradictory conditions can 
be easily reconciled, once we are reminded that mean field 
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theories do not account for discrete variables (zero not be- 
ing a 'variable' in our dynamics!). With simulations, we 
discovered complex extinction scenarios (e.g., non-trivial 
clustering in Fig. HJd) displayed by the stochastic model. 
However, to predict the properties of such distributions 
will be a serious challenge, as will be the task for comput- 
ing the essentials of P ({N m } Work is in progress to 
study a related, simpler problem: How does the distribu- 
tions of R's and Q's evolve? 

Clearly, the scope of our study is quite limited, so that 
many natural questions can be raised. Strictly cyclic com- 
petition in multiple species is rare in reality. What other 
surprises can we expect if we incorporate into our four 
species some of the other twelve possible rates? Simi- 
larly, if we introduce realistic birth/death rates (for bi- 
ological species, e.g.), will the lack of N conservation pro- 
duce novel behavior? Needless to say, the possibilities for 
generalization (e.g., to M > 4 species) are limitless, even 
for system with no spatial structure. In fact, effects sim- 
ilar to those discussed here are expected for other even 
numbers of species [28| l32] . Finally, as we noted in the In- 
troduction, qualitatively new phenomena (e.g., clustering, 
pattern formation, moving fronts) tend to emerge when a 
population dynamics is placed on some underlying spatial 
structure. For example, in spatial systems coexistence of 
all four species can be maintained even for A ^ [32] . Of 
course, in nature, different species are likely to compete in 
inhomogeneous environments. Models that include such 
realistic settings [43] will certainly display a richer variety 
of properties and will, hopefully, lead to a better under- 
standing of population dynamics. 

* * * 
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